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Metropolis Monte Carlo simulation is a powerful tool for studying 
the equilibrium properties of matter. In complex condensed-phase 
systems, however, it is difficult to design Monte Carlo moves with 
high acceptance probabilities that also rapidly sample uncorrelated 
configurations. Here, we introduce a new class of moves based 
on nonequilibrium dynamics: candidate configurations are generated 
through a finite-time process in which a system is actively driven out 
of equilibrium, and accepted with criteria that preserve the equilib- 
rium distribution. The acceptance rule is similar to the Metropo- 
lis acceptance probability, but related to the nonequilibrium work 
rather than the instantaneous energy difference. Our method is 
applicable to sampling from both a single thermodynamic state or 
a mixture of thermodynamic states, and allows both coordinates 
and thermodynamic parameters to be driven in nonequilibrium pro- 
posals. While generating finite-time switching trajectories incurs 
an additional cost, driving some degrees of freedom while allowing 
others to evolve naturally can lead to large enhancements in ac- 
ceptance probabilities, greatly reducing structural correlation times. 
Using nonequilibrium driven processes vastly expands the repertoire 
of useful Monte Carlo proposals in simulations of dense solvated 
systems. 

Metropolis-Hastings | Markov chain Monte Carlo | molecular dynamics | ex- 
panded ensembles 

Abbreviations: MC, Metropolis Monte Carlo; MD, molecular dynamics; MCMC, 
Markov chain Monte Carlo; NCMC, Nonequilibrium candidate Monte Carlo 

I n this paper, we describe a new technique for construct- 
| ing efficient Markov chain Monte Carlo (MCMC) [1] moves 
that both have high acceptance rates and also allow rapid 
transit through configuration space, greatly enhancing con- 
vergence rates in simulations of dense solvated systems. The 
Metropolis Monte Carlo [2, 3] sampling procedure is general- 
ized by using nonequilibrium processes to generate candidates 
for equilibrium simulations. Within this framework, moves 
that are efficient for an isolated part of a system but lead to 
near-universal rejection in standard Monte Carlo simulations 
of dense mixtures can be converted to nonequilibrium pro- 
cesses that generate candidates with higher acceptance prob- 
abilities. In this new procedure, the acceptance criteria is 
related to the nonequilibrium work, rather than the potential 
energy difference used in traditional Monte Carlo moves. 

Since their introduction in the mid-twentieth century, 
Metropolis Monte Carlo (MC) [2, 3] and molecular dynam- 
ics (MD) [4] simulations have become favored tools for sam- 
pling from complex multidimensional distributions, such as 
configurations of microscopic physical systems in thermody- 
namic ensembles. However, these methods can produce highly 
correlated samples, leading to slow convergence of estimated 
expectations. While MD requires the use of small timesteps 
for numerical stability and to approximate sampling from the 
desired distribution, MC simulations can, in principle, make 
use of non-local moves that accelerate mixing of the Markov 
chain. Indeed, vast improvements in efficiency have been ob- 



tained by applying cleverly constructed move sets that exploit 
physical intuition about the system under study, such as clus- 
ter moves in Potts and Ising model simulations [5, 6]. 

Designing efficient moves requires striking a balance be- 
tween rapid traversal of phase space and ensuring reasonable 
acceptance probabilities. For complex heterogeneous systems 
such as solvated biomolecules, achieving this balance remains 
challenging. Typically, efficient moves exploit physical insight 
into kinetically slow processes and energetically favorable con- 
figurations. Often, the experimenter may possess physical in- 
sight about one component in the system (e.g. a biomolecule) 
that permits the design of moves that would be efficient in the 
absence of other components (e.g. solvent), but encounter en- 
ergetically unfavorable interactions in their presence, reduc- 
ing acceptance rates to levels where standard MC provides 
no benefit. As an illustrative example, consider a bistable 
dimer — a pair of particles interacting with a potential with 
minima in compact or extended configurations, separated by 
a high barrier (see Fig. 1). For simulations of this system in a 
vacuum, a simple and effective standard MC move is to instan- 
taneously increase the interparticle distance from a compact 
to extended configuration (or conversely, to decrease the dis- 
tance from an extended to compact configuration) . When the 
dimer is immersed in a dense solvent, however, this move is 
met with near- universal rejection because solvent molecules 
overlap with proposed configurations. 

One approach that can allow unperturbed degrees of free- 
dom to relax, and hence maintain a reasonable acceptance 
rate, is to use a nonequilibrium process to generate candi- 
date configurations. Using the appropriate acceptance crite- 
rion for the final configuration will preserve the equilibrium 
distribution. In the case of the bistable dimer immersed in 
dense solvent, the extension (contraction) may be carried out 
over a finite number of increments interleaved with standard 
Metropolis Monte Carlo or molecular dynamics steps that al- 
low the solvent to reorganize to avoid overlap with the dimer 
particles. 

The basic idea of using nonequilibrium driven processes as 
Monte Carlo moves has precedents in both the statistical [7, 8] 
and chemical [9, 10, 11, 12] literature. Among the latter, 
Athenes developed "work-bias Monte Carlo" to enhance ac- 
ceptance rates in grand canonical Monte Carlo simulations [9] , 
Stern presented a scheme to sample an equilibrium mixture 
of protonation states at constant pH in explicit solvent [11] 
(though an inexact variant was proposed earlier [10]), and 
Nilmeier et al. [12] proposed the driving of a subset of degrees 
of freedom to enhance acceptance rates (using an approximate 
acceptance criteria) . Nonequilibrium processes have also been 
used to generate configurations for parallel tempering simula- 
tions [13, 14, 15]. 

Here, we unify these ideas and significantly extend the 
application of nonequilibrium moves in physical simulations. 
We present a theoretical framework, nonequilibrium candi- 
date Monte Carlo (NCMC), that is applicable to both single 
thermodynamic states (e.g. NVT, NpT, /xVT ensembles) as 
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Fig. 1. Bistable dimer potential and instantaneous MC 
moves in WCA fluid. An extension move increases the dimer exten- 
sion by Ar = +ro, while a compaction move decreases the dimer extension by 
Ar = — ro- Both move types meet with near-universal rejection when implemented 
as instantaneous MC moves in a dense WCA fluid. Note that the lower panel is only 
a cartoon — the actual described simulation is of a 3D system. 



well as mixtures of thermodynamic states (e.g. expanded en- 
semble [16, 17] simulations). Nonequilibrium proposals may 
drive a subset of degrees of freedom, the thermodynamic pa- 
rameters characterizing the equilibrium distribution, or both, 
significantly expanding the repertoire of Monte Carlo moves 
that lead to high acceptance and efficient mixing in dense 
condensed-phase systems. 



Equilibrium and Expanded Thermodynamic Ensembles 

For physical systems in equilibrium, the probability of observ- 
ing a microstate is given by the Boltzmann distribution, 



7T\(x) 



-L 



dx e 



-u x (x) 



[i] 



where x £ T denotes a microstate of the system (which may 
include coordinates, momenta, and other dynamical variables, 
such as simulation box dimensions), A denotes a set of ther- 
modynamic parameters whose values define a thermodynamic 
state, and Z\ is a normalizing constant known as the partition 
function. 

The reduced potential u\(x) depends on the thermody- 
namic ensemble under consideration [18]. For instance, in 
an isothermal-isobaric (NpT) ensemble, the reduced poten- 
tial will assume the form, 



u\{x) 



p[H(x)+pV(x) 



[2] 



which depends on the Hamiltonian H(x) (which may include 
an external biasing potential, and is presumed to be invariant 
under momentum inversion) and the system volume V(x). In 
this ensemble, the vector of controllable thermodynamic pa- 
rameters A = {/3,H,p} includes the inverse temperature f3, 
the Hamiltonian H(x), and external pressure p. Other ther- 
modynamic parameters and their conjugate variables can be 
included or excluded to generate alternative physical (or un- 
physical) ensembles. 

To allow sampling from multiple thermodynamic states 
within a single simulation, we also define an expanded ensem- 
ble [16, 17], which specifies a joint distribution for (x, A) in a 



weighted mixture of thermodynamic states, 

Z\7T\(x)uJ\ 



7r(x, A) 



E J T dyZ u 7r u (y)c 



[3] 



where uj\ > specifies an externally-imposed weight for state 
A. Here, A £ Q may assume values in a discrete or continuous 
space Q. If the set Q consists of a single value of A, a single 
thermodynamic state is sampled, and tt(x, A) = tt\(x). These 
thermodynamic states may correspond to a variety of different 
states of interest, such as temperatures in a simulated tem- 
pering simulation [19], alchemical states in a simulated scaling 
simulation [20], or protonation states in a constant-pH simu- 
lation [21]. 



Nonequilibrium Candidate Monte Carlo 

We first describe the general form of NCMC. At the start of an 
iteration, the current sample in the Markov chain, (x^ n \ A^), 
which is assumed to be drawn from tt(x, A), is used to initialize 
a trajectory, (xo, Ao) = {x^ n \ A^). A candidate configuration 
(xt,At) is then proposed through a nonequilibrium process 
in which a set of degrees of freedom and/or thermodynamic 
parameters may be driven according to some protocol [22] se- 
lected with a probability dependent only on (xo,Ao). Even 
if we only wish to sample from a single thermodynamic state 
A, we may use a protocol that transiently drives the ther- 
modynamic parameters away from A and back again (as in 
Ref. [14]). Finally, an acceptance probability is computed 
and used to decide whether the next sample in the Markov 
chain, (x^ n+1 \ \^ n+1 ^), is the candidate, (xt, At), or the mo- 
mentum reversal of the initial sample, (x^ n \ A^). 

An NCMC move begins by selecting a protocol A from 
a set of possible protocols with probability P(A|xo, Ao), such 
that there exists a reverse protocol labeled as A (to be defined 
momentarily) with P(A|xt,At) > 0. A protocol A specifies 
both a series of T perturbation kernels a t (x,y) and propa- 
gation kernels Kt(x,y), arranged in an alternating pattern 
such that A = {ai, Ki, ct2, K2, . . . , ckt, Kt}. Both a t (x,y) 
and Kt(x,y) are conditional probabilities of y £ T given any 
x £ T, and must satisfy the requirement that if p(x,y) > 0, 
then p(y,x) > 0, for p substituted by at and Kt. 

Each perturbation kernel a t drives some or all of the de- 
grees of freedom x in a stochastic or deterministic way (e.g. by 
driving a torsion angle, a distance between two atoms, or 
the volume of the simulation cell). Similarly, each propaga- 
tion kernel K t propagates some or all of the coordinates of 
the system at fixed At according to some form of MCMC or 
MD (e.g. Metropolis Monte Carlo [2, 3], velocity Verlet [23] 
deterministic dynamics, or overdamped Langevin stochastic 
dynamics [24, 25]) that may also depend on the time index 
t. Interleaving perturbation and propagation allows for ener- 
getically unfavorable interactions introduced by perturbation 
to be relaxed during propagation, potentially increasing ac- 
ceptance rates relative to the instantaneous perturbations of 
standard Metropolis Monte Carlo. 

The procedure by which a trajectory X = (xq,xi, . . . , xt) 
is generated from an initial microstate xo according to a pro- 
tocol A can be illustrated by the scheme, 



x 



X 1 > X\ 



Xt-1 > X T > XT 



[4] 



1 The described pathwise detailed balance condition is closely related to "super-detailed balance" 
(see, e.g. [26]), but also accounts for momentum reversal to extend its definition to include molec- 
ular dynamics integrators. 
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Application of the perturbation a t to xt-i generates a per- 
turbed configuration x$ , which is then propagated by the ker- 
nel K t to obtain x±. 

The reverse protocol A = {Kt, ckt, • • • , Ko, «o} reverses 
the order in which the perturbation and propagation steps 
are applied, generating the time-reversed trajectory X = 
{xt, . . • , xo}, where x denotes x with inverted momenta, 



XT > X T 



XT-1 



• — >> xi — > x 1 — > x [5J 

The next step in NCMC is to accept or reject (xt, At) as 
the next sample in the chain. To ensure that the stationary 
distribution tt(x, A) is preserved, we enforce a strict pathwise 
form of detailed balance, 1 

A(X\A) U(X\x ,A) P(A\x , A ) tt(x , A ) 

- A(X\A) U(X\x T ,A) P(A\x T , At) tt(x t , At). [6] 

The quantity A(X|A) is the NCMC acceptance probability, 
while II(X|xo,A) and U(X\xt-, A) denote the probability of 
generating trajectory X from initial configuration xo using 
protocol A, or X from final configuration xt with protocol A, 
respectively, 



U(X\x ) 
U(X\x T ) 



— ] [ a t (x t - 

l<t<T 



i Kt(xl,x t ) 
i K t (x t ,Xt) 



[7] 



After generating (xt, At) and evaluating A(X\A), we generate 
a uniform random variate U. If U < A(X\A), then the candi- 
date becomes the next value in the chain, (x^ n+1 \ X^^ 1 ^) = 
(xt,At). Otherwise, it is rejected, we perform a momentum 
flip, and the next value becomes (x^ n+1 \ \( n+1 ^>) — (x ,Ao). 
Alternately, we may perform a momentum flip upon accep- 
tance, {x^ n+1 \ A( n+1 )) = (xt,At) and preserve the momen- 
tum upon rejection, (x^ 1 ^ , A ( ^ n+1 ' ) ) — (xo,Ao). We cannot, 
however, ignore the momentum flip completely; as explained 
in the Appendix, it is necessary to preserve the equilibrium 
distribution. 

We note that NCMC need not be used exclusively to sam- 
ple from 7r(x, A), but can be mixed with other MCMC moves 
or with MD [1]. For example, one may reinitialize velocities 
from the Maxwell-Boltzmann distribution after each NCMC 
step; this is a Gibbs sampling MCMC move using the marginal 
distribution for velocities. 



Perturbation Kernels 

A large variety of choices are available for the perturbation 
kernels a t (x,y). Through judicious selection of these kernels, 
a practitioner can design nonequilibrium proposals that carry 
some component of the system from one high-probability re- 
gion to another with high acceptance rates. We briefly de- 
scribe a few possibilities. 



Summation of Eq. 6 over all trajectories starting with xo and 
ending with xt recovers the standard detailed balance condi- 
tion (see Appendix for proof). 

We define the ratio of proposal kernels as, 



a(X\A) 
a(X\A) 



at(xt,xt-i) 
tJi at(x t -i,x* t y 



n 



[91 



and the ratio of propagation kernels as the exponentiated dif- 
ference in forward and backward conditional path actions as, 



-A5(X|A) 



n 



K t {x u x* t ) 
\ K t (x*,x t )' 



[101 



Using the above expressions and the momentum invariance 
property 7r(x,A) = 7r(x,A), we may write the ratio of accep- 
tance probabilities as, 

A(X\A) _ tt(x t , At) P(A\xt, Xt) U(X\x T , A) 
A(X\A) ~ 7r(x ,Ao) P(A|x ,Ao) n(X|x ,A) 

~ t 

_ 7t(xt, At) P(A|x t , At) t-t at(x$,x t -i) K t {x t ,xX) 
7r(x ,Ao) P(A|x ,Ao) OLt(xt-\,x%) K t (xt,x t ) 



lot P(A|x T , At) a(X\A) 
W P(A|x ,Ao) a(X\A) ( 



-A5(X|A)-Au(X|A) 



[ill 



where Au(X\A) = ut(xt) — uo(xo) is the energy difference. 
Eq. 11 is the main result of this paper, and is highly general 
with regard to both the choice of protocol for perturbation 
and propagation. In subsequent sections, we discuss specific 
choices for these protocols that lead to particularly simple 
acceptance criteria. 

Many choices of acceptance probabilities A(X|A) that sat- 
isfy Eq. 11 are possible, including the well-known Metropolis- 
Hastings criterion [2, 3], 

A(X\A)= [12] 



. f uj t P(A\xt, At) a(X\A) 
mm < 1, ^Vin r nrrrrr ^ 



u P(A|x ,Ao) a(X\A) 



■A5(X|A)-An(X|A) 



Stochastically Driven Degrees of Freedom. Suppose we wish 
to drive a torsion angle <j> (an angle subtended by four bonded 
atoms) stochastically by rotating it to a new torsion angle 
cj)' (holding bond lengths and angles fixed) according to some 
probability, such as the von Mises circular distribution cen- 
tered on 6, 



vtfW = [2tt/o(/s)]- 



S(0 '-</>) 



[131 



with Io(k) denoting the modified Bessel function of order 
zero and k > taking the role of a dimensionless force con- 
stant. Because the stochastic perturbation is made in a non- 
Cartesian coordinate, a Jacobian J((j>) must be included to 
compute a(x,y) in Cartesian coordinates, resulting in the ra- 
tio, 



c*t(y,x) = y(<l>\<f>')J(<f>) 
a t (x,y) riWlMW) 



[14] 



where J{<j)') = J{<j>) — 1 because the transformation (a rota- 
tion about a bond vector) preserves the Cartesian phase space 
volume. 



Deterministically Driven Degrees of Freedom. Instead of per- 
turbing the torsion angle stochastically, we can deterministi- 
cally drive it in small, fixed increments A(j>. In this case, we 
effectively define an invert ible map M that takes x —> y, such 
that y = Mx differs from x only in the rotation of the spec- 
ified torsion by Acj). To implement this, we may choose a 
perturbation Acj) from a distribution where zLAcp have equal 
probability, and drive <j>{x) from its current value cfio to a fi- 
nal value cj)T — 0o + Acj) over T steps in equal increments, 
such that <j>{xt) is constrained to <j>t = (l — t/T)<fio + (t/T)<pT- 
In this case, ott(x,y) = S(y — A4x)J(x), where the Jacobian 
J(x) represents the factor by which Cartesian phase space is 
compressed on the application of the map which is again 
unity for rotation about a torsion angle by A<p, and, due to 
the invertibility of the map, the ratio at(y,x)/at(x,y) = 1. 
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Simulation Box Scaling. Another possible deterministic per- 
turbation kernel is simulation box scaling. A barostat can 
be implemented by proposing propagation kernels that scale 
the molecular centers and box geometry by a factor s = 
[(V(x) + AV)/V(x)] 1/s with AV chosen uniformly from [V - 
AVb, V + AVb] applied as a factor of s 1//T over the course of T 
steps. In this case, the perturbation kernel a t (x,y) is a delta 
function that compresses or expands the molecular centers 
and box geometry. Since the Jacobian is the ratio of infinites- 
imal volumes upon scaling, the ratio of perturbation kernels 
is a (X | A) /a (X | A) = s 3N , where N denotes the number of 
molecular centers. 

Thermodynamic Perturbation. In many driven nonequilib- 
rium processes, there is no direct perturbation to the co- 
ordinates, such that a t (x,y) = 5(x — y) and the ratio 
a(X\A) I f a(X\A) = 1. Instead, only the thermodynamic pa- 
rameters A are varied in time, carrying the system out of equi- 
librium through action of the K t propagation kernels. We 
recover Neal's method [7] if the reduced potential ut is a sim- 
ple linear interpolation such that u t (x) = (1 — t/T)uo(x) + 
(t/T)ur(x), the probability of choosing protocol A is symmet- 
ric with A, and MC [2, 3] is used for the propagation kernel 
K t . 
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Propagation Kernels 

The choice of propagation kernels available is also very broad. 
If strong driving is performed in selection of a, one may elect 
to choose a time-independent propagation kernel K t (x,y) = 
K(x,y) that samples from a stationary distribution tv(x) of 
interest. Alternatively, a strongly time-dependent K t could 
be selected to transiently drive the system out of equilibrium, 
or from the equilibrium distribution at one thermodynamic 
state to another. Some possible choices are described below. 

Reversible Markov Chain Monte Carlo. One may propagate 
some or all of a system's degrees of freedom (e.g. those not 
affected by the perturbation kernel at) by a method that sat- 
isfies detailed balance in 7r*, 

7Y t (x)K t (x,y) = 7T t (y)K t (y,x), [15] 

where 7r t (x) = Z^ 1 e~ Ut ^ for a specified u t (x). Many 
MCMC methods [1], including Metropolis [2, 3] and various 
hybrid Monte Carlo (HMC) algorithms that combine discrete- 
timestep integrators with Monte Carlo acceptance/rejection 
steps [27, 28], obey detailed balance and are commonly used 
for the simulation of physical systems. 

By analogy with Crooks [29] , we define a work w and heat 
q for the nonequilibrium driven process, 

T 

w{X\K) = ^2[ut(xt)-ut-i(x t -i)] [16] 
t=i 

T 

q{X\K) = -ut(s t *)] [IT] 

t = l 

such that w(X\X) + q(X\X) = Au(X\X), a restatement of the 
first law of thermodynamics. 

The conditional path action difference can then be written 
in terms of the heat of the process, q(X\A), 

AS(X\A) = -Infl^p- = -q(X\A), [18] 



Fig. 2. Trajectories of WCA dimer system in vacuum and 
solvent. Left: The dimer extension r as a function of simulation iteration. The 
dotted horizontal line denotes division between compact and extended configurations. 
The quantity r printed above each plot indicates the estimated integrated autocorre- 
lation time for the dimer extension r. Right: Histogram accumulated over trajectory 
(black), with true equilibrium distribution (red). Plot titles denote whether simulation 
was run in vacuum [vacuum) or dense WCA fluid [solvent), and whether the sim- 
ulation utilized only 500 GHMC steps per iteration [MD) or included instantaneous 
MC [MC) or 2048-step NCMC moves [NCMC) following each iteration. 

leading to an acceptance probability similar to standard MC, 
except that the work, w(X\A), replaces the instantaneous po- 
tential energy difference, 

A(X\A) = uj t P(A\x T ,\T)a(X\A) c _ w(xlA) 
A(X\A) ujo P(A|x ,Ao) a(X\A) 6 ' [ 1 

Deterministic Dynamics. When an isolated system is propa- 
gated by a symplectic integrator — a reversible, deterministic 
integrator that preserves phase space volume — the propaga- 
tion kernels follow K t (x,y) = K t (y,x). Hence, AS(X\A) = 
and the acceptance ratio is, 

A(X\A) = ujt P(A\xt,\t) a(X\A) Au(xlA) 
A(X\A) cjo P(A|xo,A ) a(X\A) ' L J 

where Au(X\A) = ut(xt) — uo(xo) is the energy difference. 
The equivalence of the work and energy difference for volume- 
preserving integrators was previously recognized in the con- 
text of fluctuation theorem calculations [30, 31]. 

Symplectic integrators include velocity Verlet [23] . These 
integrators are also symplectic when utilizing constraints 
through the application of algorithms such as RATTLE [32], 
provided that the constraints are iterated to convergence each 
timestep [33]. 



Alternatives to using NCMC to correct stochastic integration include introducing a Metropoliza- 
tion correction after each step, as in the generalized hybrid Monte Carlo (GHMC) integrator we use 
in the example [37, 1, 36, 28]. 
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Fig. 3. Acceptance probabilities of NCMC proposals. Top: 
Acceptance probability of NCMC proposals as a function of length of nonequilibrium 
proposal trajectory (black dots), compared with instantaneous MC proposal (red line). 
Inset: Enlarged region with acceptance probability shown on linear scale. Estimated 
95% confidence intervals are shown as vertical lines. 



16 32 64 128 256 512 1024 2048 4096 8192 
NCMC switching steps 




Stochastic Dynamics. Stochastic integrators such as velocity 
Verlet discretizations of Langevin dynamics [34, 35] sample 
a modified distribution that differs from the desired equilib- 
rium distribution 7r* in a timestep-dependent manner [36]. 
While this modified distribution may be difficult or im- 
possible to compute in order to recover equilibrium prop- 
erties by reweighting, computation of the relative action 
A«S(X|A) is relatively straightforward, and the NCMC ac- 
ceptance criteria ensures that the NCMC-sampled configu- 
rations are distributed according to the desired equilibrium 
ensemble 2 . As examples, we compute A<S(X|A) for the over- 
damped Langevin (Brownian) dynamics integrator of Ermak 
and Yeh [24, 25] and the Briinger-Brooks-Karplus (BBK) 
Langevin integrator [38, 39, 40] in the Appendix. 

Illustrative Application: Bistable Dirtier in a WCA Fluid 

To demonstrate NCMC, we ran simulations of a bistable 
dimer (adapted from Section 1.3.2.4 of Ref. [36]) in vacuum as 
well as a dense fluid. The dimer consists of a pair of "bonded" 
particles interacting via a double-well potential, with minima 
at r = ro (compact) and r = 2ro (extended), and a 5/cbT 
barrier (see Fig. 1). In the solvated simulations, the dimer 
was immersed in a dense bath (reduced density pa 3 = 0.96) 
of particles that interact with the bonded particles and each 
other via the Weeks-Chandler- Andersen (WCA) soft repulsive 
potential [41]. Each simulation "iteration" consisted of veloc- 
ity reassignment from the Maxwell-Boltzmann distribution, 
500 steps of generalized hybrid Monte Carlo (GHMC) dynam- 
ics [37, 1, 36, 28] (essentially, a Metropolis-corrected form of 
Langevin dynamics, henceforth referred to here as MD), op- 
tionally followed by either an instantaneous MC move or an 
NCMC move. 

The rate at which effectively uncorrelated samples are gen- 
erated can be quantified in terms of the correlation time r for 
the dimer extension r{t) (shown in Fig. 2). This time repre- 
sents the asymptotic decay time constant for the correlation 
function C(t) = (r(O)r(t)), which will behave like 

C(t) = Coo + (Co - Coo)e~ t/T [21] 

for large t, where Co = (r 2 ) and Coo — (r) 2 . The correla- 
tion time is related to the statistical inefficiency, g — 1 + 2r, 
a factor that describes the number of iterations necessary to 
generate an effectively uncorrelated sample [42]. 

For the MD simulation in vacuum (Fig. 2, top trace), we 
observe slow hopping between compact and extended config- 
urations with a correlation time r = 59.2 iterations, resulting 
in slow convergence of the histogram. Introducing instanta- 
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NCMC switching steps 

Fig. 4. Statistical efficiency gain of NCMC proposals rel- 
ative to instantaneous MC proposals. Top: Effective correlation 
time r e ff , in iterations, for MD+NCMC (black dots) compared to MD alone (red 
line). Bottom: Relative statistical efficiency of MD+NCMC, in terms of number of 
uncorrelated configurations generated for a fixed amount of computational effort, for 
MD+NCMC (black dots) relative to MD alone (red line). 

neous MC dimer extension/contraction moves that modify the 
dimer extension by Ar = =bro reduces the correlation time to 
r ~ 0.0, such that an uncorrelated configuration is generated 
after each iteration of 500 MD steps and one instantaneous 
MC step (Fig. 2, second trace from top). 

When the dimer is immersed in a dense fluid of WCA par- 
ticles, however, iterations consisting of 500 MD steps alone 
result in extremely slow barrier crossings, requiring g « 600 
iterations to produce an uncorrelated sample (Fig. 2, middle 
trace). Unlike in vacuum, the introduction of instantaneous 
MC moves does not significantly reduce the correlation time 
r (Fig. 2, second trace from bottom). However, perform- 
ing these same dimer expansion and contraction moves over 
2048-step NCMC moves (Fig. 2, bottom trace) allows the 
system to rapidly mix between both compact and extended 
states with a correlation time of r = 4.0 iterations. While 
each iteration requires a 5-fold increase in computational ef- 
fort (500 MD steps + 2048 NCMC switching steps = 2548 
force evaluations, versus 500 force evaluations for MD alone), 
a 67-fold reduction in correlation time is achieved, resulting 
in a remarkable order-of-magnitude gain in overall efficiency. 

The length of the NCMC switching process is a free pa- 
rameter that may be tuned to further improve efficiency. To- 
wards this end, we estimated the acceptance probability of 
the extension/contraction moves in dense solvent as a func- 
tion of switching length (Fig. 3). While instantaneous MC 
proposals of =bro are only accepted with probability « 10 -27 
(the error is this quantity is likely underestimated due to its 
extremity), dividing the move into smaller steps boosts the 
acceptance rate to a level useful in condensed-phase simula- 
tion. If the move is divided into a small number of steps (1 to 
8), there is little to no increase in acceptance rate, but for an 
intermediate number of steps (16 to 1024), there is a superlin- 
ear boost in the acceptance probability relative to the length 
of the switching process. The acceptance probability finally 
reaches useful levels around 2000 steps, achieving an accep- 
tance rate of 12% using nonequilibrium proposal trajectories 
of 2048 steps, or 38% for 8192 steps. 
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In general, there is no direct relationship between accep- 
tance probability and efficiency. Under certain assumptions 
relevant to the bistable dimer, however, it is possible to link 
the NCMC acceptance probability to r e ff , an indirect estimate 
of the correlation time, 



T e ff = TMD 



tncmc 



tmd + tncmc 



[221 



where tmd and tncmc are correlation times for iterations con- 
sisting solely of MD or NCMC moves, respectively. The latter 
correlation time may be estimated from the average accep- 
tance probability 7 by tncmc ~ — 1/ ln(l — 27) (see Appendix 
for derivation). 

As shown in Fig. 4, the effective correlation time r e ff is 
only diminished when the NCMC acceptance probability is 
large enough such that tncmc ~ tmd, which occurs when 
7 > 0.13% (about 256 switching steps or more). For shorter 
switching times, even though the acceptance probability is 
high relative to instantaneous MC, it is still too small to sig- 
nificantly reduce r e ff. 

When comparing efficiency, we are most interested in 
the rate of generating uncorrelated configurations for a given 
amount of computational effort. Relative to MD alone, this 
rate is described by the efficiency gain, 



#ncmc(7md + Tncmc) ' 



[231 



Here, Tmd = 500 steps per iteration, and Tncmc is again var- 
ied from 1 to 8192 steps. The results are shown in the bottom 
panel of Fig. 4. Surprisingly, there is actually a slight loss in 
efficiency at short switching times — dropping to a minimum 
of 86.9% the efficiency of MD alone at 128 steps — followed by 
a rapid gain in efficiency, plateauing at an efficiency gain of 
~ 13 X the efficiency of MD alone for 2048-4096-step NCMC 
proposals. (A similar plateau behavior is observed in the mod- 
ified parallel tempering protocol of Ref. [15] .) After this point, 
longer switching times do not achieve as high of an efficiency 
gain; even though the acceptance rate continues to increase 
as the number of NCMC switching steps is doubled again to 
8192 steps, the reduction in correlation time is not sufficient 
to offset the additional cost of these moves. 



Epilogue 

We have described a procedure — nonequilibrium candidate 
Monte Carlo (NCMC) — by which nonequilibrium proposals 
can be used within MCMC simulations to enhance accep- 
tance rates. In our illustration, we have demonstrated how its 
use can lead to large improvements in statistical efficiency — 
the rate at which uncorrelated samples are generated for a 
fixed amount of computational effort. In other applications, 
whether similarly large efficiency gains are achieved will de- 
pend on the precise nature of the system under study and the 
nonequilibrium proposals introduced. The most straightfor- 
ward approach — borrowing Metropolis Monte Carlo proposals 
that are reasonable for one component of the system in iso- 
lation, and converting them to nonequilibrium proposals — is 
likely to be a fruitful avenue for generating efficient schemes, 
as was demonstrated here for a simple system. 

More generally, the problem of selecting efficient nonequi- 
librium proposals is similar to the problem of choosing good 
reaction coordinates, in that it is desirable to drive the 
system along (possibly complex) slow collective coordinates 
where orthogonal degrees of freedom relax quickly. The 
search for such collective coordinates is a topic of active re- 
search [43, 44, 45, 46, 47, 48, 49]. Given an initial nonequilib- 



rium protocol, the issue of optimizing such a protocol to min- 
imize dissipation (and maximize acceptance) is also a topic of 
active study, led by forays into the world of single-molecule 
measurement [50, 51, 52]. Recent work has also suggested 
that estimating the thermodynamic metric tensor along the 
nonequilibrium parameter switching path [53, 54, 55, 56], 
could prove useful in adaptively optimizing the switching pro- 
tocol [57]. 

We note that switching trajectories contain potentially 
useful information. Indeed, several methods [58, 59, 56] have 
recently been developed to estimate equilibrium properties 
from nonequilibrium samples through the application of sta- 
tistical estimator theory to nonequilibrium fluctuation theo- 
rems [30, 60, 61]; these are particularly relevant to switch- 
ing between multiple thermodynamic states. Though the de- 
velopment of efficient estimators that utilize both nonequi- 
librium switching trials and sampled equilibrium data gener- 
ated in NCMC simulations remains an open challenge, it is 
at least straightforward to incorporate information from re- 
jected NCMC proposals in the estimation of equilibrium av- 
erages [26, 62]. 

Materials and Methods 



WCA Dimer Simulations. The dimer system considered here consists of two parti- 
cles that interact via a double-well "bonded" potential in the interatomic distance r, 

2 



Uhond(r) 



h 



(r - r - s) 



[24] 



Uwca(v) 



4e 



with h = 5 k&T, To — TwCA, and S — TwCa/2, where TwCA = 2 1 ^ 6 (7. 
Simulations denoted as "vacuum" contain only these two particles, while simulations 
denoted as "solvated" also interact with a dense bath of particles via the WCA non- 
bonded potential [41], 

[(7) 12 -(7) 6 ]+c, ^<nvcA )[25] 

,0 r > r W CA 

with mass m — 39.9 amu, G — 3.4 A, and 6 = 120 k^,T . The nonbonded 
WCA interaction is excluded between the two bonded particles. The solvated system 
contains a total of 216 WCA particles at a reduced density of pa 3 = 0.96. For 
all simulations, the reduced temperature is kj$T j € — 0.824. A custom Python 
code making use of the GPU-accelerated OpenMM package [63, 64, 65] and the Py- 
OpenMM Python wrapper [66] was used to conduct the simulations. All scripts are 
available for download from http : / / simtk . org/home/ncmc. 

To ensure that observed differences were not due to changes in the sta- 
tionary distribution of the integrator, we used generalized hybrid Monte Carlo 
(GHMC) ([37, 1, 36, 28]) for all our simulations. GHMC is based on a velocity 
Verlet discretization [23] of Langevin dynamics — the two are equivalent in the limit 
of small timesteps - but includes an acceptance/rejection step to correct for errors 
introduced by finite timesteps so that the stationary distribution is the exact equilib- 
rium distribution. We used a timestep of 0.002r, where T — yj a 2 m/e, and 
the collision rate was set to T~ 1 . With this timestep, the acceptance probability is 
99.929i0.001%; the resulting dynamics closely approximates Langevin dynamics. 

In simulations employing instantaneous Monte Carlo moves, a perturbation Ar 
to the interatomic distance V of the dimer was chosen according to, 



Ar 



if r < 1.5ro 

if 1.5ro < r < 3ro . 

otherwise 



[26] 



The dimer was contracted or expanded about the bond midpoint to generate a new 
configuration X new with dimer extension r n ew from the old configuration X Q ld 
with dimer extension r \d, and the move accepted or rejected with the Metropolis- 
Hastings criterion [3], 

-0[tf(aJnew) + tf(*old)] T (<r u v 



Fold) 



-)} 
[27] 
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where the Jacobian ratio term J r (x ld, ^new) — (^new/^old) 2 accounts for 
the expansion and contraction of phase space due to the Monte Carlo proposals. 

For simulations employing T-step NCMC moves, proposals were made by se- 
lecting a new velocity vector from the Maxwell-Boltzmann distribution, integrating T 
steps of velocity Verlet dynamics [23] for all bath atoms as the dimer extension was 
driven from r Q ld to r ne w in equal steps of size Ar/T, and accepting or rejecting 
based on the modified Metropolis criteria for symplectic integrators (Eqs. 12 and 20), 



A(X) = min{l,( 



-/3[H(x T )-H(x )} 



J r (x ,x T )}. [28] 



The Jacobian ratio is also (r n ew/^old) 2 ■ This scheme corresponds to a choice for 
the perturbation kernel of, 



att(x,y) = 



r(y) 



r(x) 



6([r(y)-r(x)]-[Ar/T}), 



[29] 



where r(x) denotes the dimer separation of configuration X. The propagation kernel 
Kf(x, y) corresponds to velocity Verlet dynamics where the dimer atoms are held 
fixed in space. The final configuration after the MC or NCMC rejection procedure 
was stored and plotted to generate Fig. 2. 

The mean acceptance probabilities for each switching time T can be estimated 
via the sample mean 



1 N 



N 



[30] 



For numerical stability, logarithms of A(X n ) were stored, as CL n = In A(X n ) 
We then estimated In (A) (shown in Fig. 4) as 



\n{A) T n -lniV + lnfr+^V" 

n=l 

vhere b = max n CL n . 



[31] 



Integrated autocorrelation times were estimated using the rapid scheme described 
in Section 5.2 of Ref. [42]. 

The acceptance probabilities plotted in Fig. 4 were estimated from a trajectory 
consisting of 10 000 iterations of 2048-step NCMC, with 500 steps of GHMC dy- 
namics in between each NCMC trial, to ensure equilibrium sampling. Prior to each 
2048-step NCMC iteration, a velocity from the Maxwell-Boltzmann distribution was 
selected, and NCMC trial moves with varying switching times were conducted solely 
to accumulate statistics. The statistical error in the estimate of In (A) T and the 
computed relative efficiency over instantaneous Monte Carlo was estimated by 1000 
bootstrap trials, in which the dataset of 10 000 work samples was resampled with 
replacement in each bootstrap trial and 95% confidence intervals computed from the 
distribution over bootstrap replicates. 

The reference distribution for the interparticle distribution 'P(r) plotted in red 
on the right side of Fig. 2 was computed analytically for the vacuum system, 



^vac(r) 



47rr 2 e 



PU hond (r) 



[32] 



For the solvated system, this distribution was estimated from an umbrella sampling 
simulation employing a modified bonded potential intended to remove the barrier in 



between compact and extended states, 

^umbreiia(r) = k B T\nr 2 + #(r min - r)(K/2)[r - r min ] 2 
+ 0(r - r max )(K/2)[r - r max ] 2 , [33] 



where r m i n — ro, r max — 2.05ro, and K — k B T ' I 'rj 2 , with r\ — 0.3 A, 
and 0(r) is the Heaviside function that assumes a value of unity for r > and 
zero otherwise. The true solvated interparticle distribution p(r) was estimated by 
reweighting the data produced from this simulation, using the relationship 



Psoi(r) 



~(3 [C/bond (r n ) - *7 umbr ella On)] 



^umbrella 



[34] 



where r n denotes the bond separation for sample n, and a finite-width histogram 
bin was used instead of the delta function S(r). 
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Fig. 5. Umbrella sampling simulation of the dimer in WCA 
solvent. Left: The dimer extension r as a function of simulation iteration. 
Right: The histogram accumulated over the trajectory with the observed histogram 
in black and the reweighted histogram (corrected for the applied umbrella bias poten- 
tial) in red. 
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Appendix 



Proof that NCMC preserves the equilibrium distribution 

Following the proof for GHMC in Ref. [36], here we show that NCMC preserves the equilibrium distribution. The expected 
acceptance rate for NCMC moves initiated from (x, A) is, 



a(x,X) = J dA J dX p(X,A\x,X)A(X\A). [35] 



Suppose that we have a variate (x^ n \ A^) drawn from the equilibrium distribution ir(x, A). The probability density of the 
next value in the chain, p(x < " n+1 \ A ( - n+1 ' ) ), has contributions from two scenarios: when the candidate variate is rejected and 
when it is accepted. The contribution from rejecting the candidate and flipping the momentum such that (x^ n+1 \ X^ n+1 ^) = 
(x (n) ,A (n) ) is, 



[361 



f dx < x > A )[! - <*(x, X)]S(x - x {n+1) )6(X - A (n+1) ) = 7r(x (n+1) , A (n+1) )[l - a(x (n+1) , A (n+1) )]. 
J x 

The latter contribution from accepting the candidate such that (x ( ^ n+1 ' ) , A ( ^ n+1 ' ) ) = (it, At) is, 

J dx^V(x,A) J dX J dAp(X,A\x,X)A(X\A)S(xT-x (n+1) )6(X T -X (n+1) ) 

= J dxoj^ J dx J dA [7t(x ,Xo)p(X,A\xo,Xo)A(X\A)]6(xt - x (n+1) )6(X T - A (n+1) ) 

- f+rprtfa Art MX A, iT , A T , MM)} H*r - rt"«>) 4 (A. - A<->, 

= ir(x (n+1 \\ (n+1) )a(x (n+1 \\ (n+1) ), [37] 

where p(X, A\xo, Ao) = II(X|xo, A)P(A|xo, Ao) is the probability of generating the trajectory-protocol pair (X, A) from (xo, Ao), 
and the pathwise detailed balance condition (Eq. 6) is used to produce the quantity in brackets. 
Taking the sum of Eqs. 36 and 37, we find that the equilibrium distribution is preserved, 

p(x (n+1) ,A (n+1) ) =7r(x (n+1) ,A (n+1) ) [38] 

By analogous reasoning, maintaining the momentum upon rejection, {x^ n+1 \ X^ n+1 ^ — (x^ n \ A^), and flipping it upon 
acceptance, {x^ n+1 \ X^ n+1 ^ — (xt, At) will also preserve the equilibrium distribution. 



Acceptance criteria for overdamped Langevin (Brownian) integrator of Ermak and Yeh 

A common integrator for Brownian dynamics (the overdamped regime of Langevin dynamics) , in which only coordinates x are 
explicitly integrated, is given by Ermak and Yeh [24, 25]. In our notation, where the perturbed coordinates Xt are propagated 
by one step of the stochastic integrator to obtain xt, application of the propagation kernel K(xt,x t ) can be written, 



x t = 



At 



7m \77Ti ) 



[391 



where m is the particle mass, F t (x) = —(d/dx)H t (x) is the (potentially time-dependent) systematic force, and 7 is an effective 
collision frequency or friction coefficient with units of inverse time. The noise history £t for each degree of freedom is a normal 
random variate with zero mean and variance drawn from the distribution 



exp 



2 4t 



[401 



In NCMC, every application of the propagation kernel Kt produces a transition x^ — > xt determined by the noise history 
variable there is a corresponding £t that generates the opposite step, xt —^x*. By noting 



we can compute the relationship, 



Footline Author 



X t = 



X t 



X t 



—F t (x* t ) + V2(—Y & 
\imj 

{ A/\ 1/2 
F t (x t ) + V2( ^ ) ^ 



7m 

At 

7m 



\7m J 



[41] 



1 /AA 1/2 

= 71 ^ [F t ( Xt ) + F t ( X ;)]-t t , [42] 
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Then, the ratio of transition kernels appearing in Eq. 10 can be written in terms of noise history £t and the computed 
reverse noise history 

where the tildes are dropped because the microstate x contains no momenta. The quantity \dxt/d£t\ represents the Jacobian 
for the change of variables from the to xt, and the Jacobians in the numerator and denominator cancel. The quantity in Eq. 
43 can easily be computed during integration. 



Acceptance criteria for Langevin integrator of Brooks, Briinger, and Karplus (BBK) 

The Briinger-Brooks-Karplus (BBK) stochastic integrator [38, 39] is a popular integrator for simulating Langevin dynamics. 
In our notation, where the perturbed coordinates Xt are propagated by one step of the stochastic integrator to obtain xt, 
application of the propagation kernel K(x*,xt) can be written, 



v t 



r t = 



* , At 

rt+Atv't 
1 



27m 
~Ar 



1 + 



jAt 
2 



Vt + 



At 
2m 



F t (r t ) + 



27m 
~At 



[441 



where we have used a velocity Verlet discretization of the BBK integrator [40, 36]. Here rt and Vt denote the respective 
Cartesian position and velocity components of the microstate xt, 7 the effective collision frequency with units of inverse time, 
and m the particle mass. v' t is an auxiliary variable used only in simplifying the mathematical representation of the integration 
scheme. Note that we require two random variates, £ t and per degree of freedom per timestep in order for this scheme to 
be able to generate both the forward trajectory X and its time-reverse X (see, e.g., Section 2.2.3.2 of [36]). 

The noise history terms £t and £' t are normal random variates with zero mean and variance Their joint distribution 

can therefore be written, 



m^t) 



T exp 



[45] 



For every step Xt — >• x t , the positions and velocities undergo a transition (r% , v* ) — >> (rt,v t ) determined by the noise variables 
(Ct j CO- A corresponding choice of noise variables (it, it) will generate the reverse step, x t —> x*, carrying (r t , —v t ) —> (r*, —v*) 
With a little algebra, it is seen that, 



it — £t - ^/2jmAt v t 
it €t - yj2^mAtVt- 



[461 



In order to write the ratio of transition kernels appearing in Eq. 10 in terms of noise variables (£t, £t) and the computed reverse 
noise variables (£t,£t)> we must first compute the Jacobian J(£t,£t) because the random variates are not in Cartesian space, 



det 



dr t 
or t 



ovt 



which can be shown to be independent of and £' t . The conditional path action difference can now be computed, 



AS(X) 



inn 



K t (x t ,xt) 



mn 



m,i't)j(it,i' t ) 



Ep 2 + r)-(6 2 + ^)] 



11 K t (x* t ,x t ) -i = i m,Ct)j(^c t ) 

where the ratio of Jacobians J (it, it)/J(£t, Ct) cancels because the Jacobians are independent of the noise variates. 



[47] 



[481 



Derivation of effective correlation time for mixed MD/NCMC sampling 

For simplicity, we make the assumption that the system of interest has two long-lived conformational states of equal population 
with dimer extensions r c and r e . This assumption holds to good approximation for the WCA dimer example considered here, 
and may apply to other systems of interest as well. We assume that the correlation time for a fixed number of MD simulation 
steps per iteration is given by tmd, and describe the probability of finding the system ends up in a given conformational state 
after one iteration by a 2 x 2 column-stochastic transition matrix Tmd, 



1-a 

a 



a 
1 — a 



[491 
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For a 2 x 2 system whose time evolution is governed by the column stochastic transition matrix T, we can write the 
autocorrelation function for the dimer extension r as 



C(nAt) = (r(O)r(nAt)) 

= [r c r e ]T n 

= [r c r e ] U 

— (Co — Cog)^ 7 ' + Co 



"1/2 




r c 


1/2 




r e 



' 1 


u- 1 


■ 1/2 







r c 


n n 





1/2 _ 




r e 



[501 



where the transition matrix T has unitary eigenvalue decomposition UMU -1 , and \i is the non-unit eigenvalue of T. The 
constants are C = (l/2)(r% + r\) and Coo = (1/4) (r c + r e ) 2 . 

Relating this to the autocorrelation time r estimated from a timeseries, intended to reflect the fit to 



C(t) = (Co - Coo)e- 



-t/T 



[511 



we can see that r = — 1/ln/x. We then determine that the correlation time tmd = — 1/ln/iMD, with /xmd = 1 — 2a. 

Similarly, we can write the probability that the NCMC step will carry the system from one conformational state to another 
in terms of the acceptance probability 7, which we assume to be symmetric, 



Tncmc 



1-7 

7 



7 

1-7 



[521 



where we have correlation time tncmc — — 1/ln/iNCMC and /xncmc = 1 — 27. 

The effective transition matrix T e ff for iterations consisting of MD simulation steps followed by an NCMC trial is given by 



Teff = T 



MDlNCMC — 



1 — a 

a 



a 
1-a 



1-7 

7 



7 
1-7 



1 — (a + 7) (a + 7) — 2aj 
(a + 7) — 2cry 1 — (a + 7) 



where the effective correlation time r e ff 
and 7 = (1 - e" 1/TNCMC )/2, we obtain 



-1/ln/ieff, with /i e ff = 1 — 2 [(a + 7) — 2aj]. Substituting in a 



[53] 



T e ff 



In [1 — (1 — e _1 / T MD) — (1 — e _1 / r NCMc) + (1 — e -1 / r MD)(l — e _1 / r NCMC) 
1 1 tmd tncmc 



In \e~ 1 / T MDe — 1 / t ncmc1 -t- - 1 _i_ -r - 1 

111 L e e J 'MD ^ T NCI\ 



tmd + tncmc 



[541 



As a check of the accuracy of Eq. 54, we note that for MD with 2048-step NCMC switching, we compute r e ff ~ 4.0 iterations, 
using only tmd = 299.8 iterations and the NCMC acceptance probability 7 = 12.1%. The actual correlation time measured 
from a 10 000 iteration simulation is computed to be r e ff = 4.0. 
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